function [EC] = APEserdsRecon(datasetname, setoff, setofftype, setplot, vargin)

% [EC] = APESERDSRECON(datasetname, setoff, setofftype) 
%
% Reconstructs Raman spectrum from SERDS difference spectrum
%
% INPUT     datasetname :   filename to process (must contain SERDS data
%           setoff      :   setoff between both Laser lines
%           setofftype  :   0=[setoff in pixel], 1=[setoff in wavelength]
%           setplot     :   set to 0 to disable plot
%
% OUTPUT    EC          :   error code (0 = success, 1 = error)
%

%  (c) F. Dworkowski, 2012
%      Version 1.0, 2012-07-09

% ----- Check for input and adjust if empty -----
if nargin < 1
    help APEserdsRecon
    EC = 1;
    return
elseif nargin < 2
    setoff = 0;
    setofftype = 0;
    setplot = 1;
elseif nargin < 3
    setofftype = 0;
    setplot = 1;
elseif nargin < 4
    setplot = 1;
end

% ----- Try to load dataset -----
try dataset = load(datasetname);
catch ME
    error(ME.message)
    EC = 1;
    return
end

% ----- Check for existanece of SERDS data -----
if ~isfield(dataset,'SERDSdata')
    EC=1;
    error('Data does not contain SERDS data!');
    return
end

% ----- Override input with file values -----

if setoff == 0
    if isfield(dataset.SERDSdata, 'setoff')
        setoff = dataset.SERDSdata.setoff;
        if isfield(dataset.SERDSdata, 'setofftype')
            if strcmp(dataset.SERDSdata.setofftype,'Wavelength')
                setofftype = 1;
            else
                setofftype = 0;
            end
        end
    else
        EC = 1;
        error('setoff can not be 0 and was not found in file for correction')
        return
    end
end

%% ----- Compute -----
S = zeros(length(dataset.SERDSdata.datay),1);
switch setofftype
    case {'Wavelength', 1}
        lowoff = min(dataset.SERDSdata.datax)+setoff;
        lowidx = find(dataset.SERDSdata.datax>lowoff,1,'first');
        for i=(lowidx+1):length(dataset.SERDSdata.datay)
            S(i) = S(i-lowidx) + dataset.SERDSdata.datay(i);
        end
    case {'Pixels', 0}
        for i=(setoff+1):length(dataset.SERDSdata.datay)
            S(i) = S(i-setoff) + dataset.SERDSdata.datay(i);
        end
    otherwise
        error('Type of SETOFF not defined!')
        EC = 1;
        return
end

S = S+abs(min(S));                  % Correct for negative Y values
if setplot == 1
    try 
        close(468)
    end
    figure(468)
    subplot(3,1,1)
        if ~isfield(dataset,'dataset1')
            plot(1,1)
            text(0.5,1,'No Single-Laser datasets available')
        else
            plot(dataset.dataset1.datax, dataset.dataset1.datay,'r')
            hold on
            plot(dataset.dataset2.datax, dataset.dataset2.datay,'g')
            hold off
            title('Raw spectra')
            axis tight
        end
    subplot(3,1,2)
        plot(dataset.SERDSdata.datax, dataset.SERDSdata.datay)
        title('SERDS spectrum')
        axis tight
    subplot(3,1,3)
        plot(dataset.SERDSdata.datax,S,'r')
        if isfield(dataset,'dataset2')
            hold on
            plot(dataset.dataset2.datax, dataset.dataset2.datay,'g')
            hold off
        end
        title('Reconstructed spectrum')
        axis tight
end
%% ----- SAVE data -----
dataset.SERDSdata.reconsty = S;
save(datasetname,'-struct','dataset')
disp([datasetname,' was saved with reconstructed SERDS spectrum'])

EC = 0;     % Exit w/o error